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ABSTRACT 



Context. Time-dependent cooling processes are of paramount importance in the evolution of astrophysical gaseous 
nebulae and, in particular, when radiative shocks are present. Given the recent improvements in resolution of the 
observational data, simulating these processes in a more realistic manner in magnetohydrodynamic (MHD) codes will 
provide a unique tool for model discrimination. 

Aims. The present work introduces a necessary set of tools that can be used to model radiative astrophysical flows in 
the optically-thin plasma limit. We aim to provide reliable and accurate predictions of emission line ratios and radiative 
cooling losses in astrophysical simulations of shocked flows. Moreover, we discuss numerical implementation aspects to 
ease future improvements and implementation in other MHD numerical codes. 

Methods. The most important source of radiative cooling for our plasma conditions comes from the collisionally-excited 
line radiation. We evolve a chemical network, including 29 ion species, to compute the ionization balance in non- 
equilibrium conditions. The numerical methods are implemented in the PLUTO code for astrophysical fluid dynamics 
and particular attention has been devoted to resolve accuracy and efficiency issues arising from cooling timescales 
considerably shorter than the dynamical ones. 

Results. After a series of validations and tests, typical astrophysical setups are simulated in ID and 2D, employing 
both the present cooling model and a simplified one. The influence of the cooling model on structure morphologies can 
become important, especially for emission line diagnostic purposes. 

Conclusions. The tests make us confident that the use of the presented detailed radiative cooling treatment will allow 
more accurate predictions in terms of emission line intensities and shock dynamics in various astrophysical setups. 



Key words. Radiation mechanisms: thermal - Line: formation 
- Methods: numerical - Magnetohydrodynamics (MHD) 



(ISM:) Herbig-Haro objects - ISM: jets and outflows 



1. Introduction 

Astrophysical gases emit thermal radiation while under- 
going dynamical transformations. There are cases where 
the gas is so diluted that the typical timescales for cooling 
greatly exceed the dynamical ones but, in many instances, 
cooling and the related ionization/recombination processes 
for the emitting species become comparable to, or faster 
than, the dynamical evolution of the system and should 
be considered. Classical examples of intensively radiating 
gases are Hn regions, planetary nebulae, supernova rem- 
nants and star forming regions. Thus, when studying gas 
flows in such environments, particular care must be taken 
to treat the interplay between dynamics and radiation in 
the correct way. This is particularly true and crucial when- 
ever radiative shocks are involved. When cooling is strong, 
the ionization fraction of the emitting species is far from 
that at equilibrium, but evolves so rapidly with time that 
it must be treated in a time-dependent fashion. 

Under these conditions, the magnetohydrodynamic 
(MHD) equations, coupled with the equations describing 
the evolution of the emitting species and the radiative 
losses, must be solved by numerical means. The numerical 
problems posed by radiative cooling arc particularly chal- 
lenging whenever the advance time step of the integration 
is controlled by radiation/ionization rather than dynamics. 



This may happen in few grid points in the computational 
domain where the radiative losses are intense, e.g., right 
behind a shock front, where the cooling time becomes very 
small. Therefore, it is necessary to devise strategies able to 
deal with very different integration timescales: this is one 
of the points we address. 

Time dependent ionization calculations were previously 
performed for gaseous nebulae by Marten & Szczerba 
(|1997p in hydrostatic conditions. Their approach is simi- 
lar to our implementation of the ionization state treatment, 
however the ion species and implemented physical processes 
are in part different. Also, radiative cooling in optically-thin 
plasmas was previously investigated, and synthetic cooling 
functions were designed (e.g., Schmutzler & Tscharnutter 
I1993[) . Among the radiative numerical codes employed in as- 
trophysics, one can quote the hydrocode YGUAZU (Raga 
et al. [^UnU|) . ASTROB EAR M HD code (Poludnenko et al. 
120051 Berger & LeVeaue ll998p and Virginia Hydrodynamics 
- 1 (VH-1, Sutherland et al. [20031 Blondin & Lufkin [19931) . 
The MHD simulation code we use for our astrophysical ap- 
plications - PLUTO - is a freely distributed application 
developed and maintained at the Turin University - Turin 
Astronomical Observatory (Mignone et al. I2007p . A pre- 
vious numerical analysis about the evolution of radiative 
shocks in Young Stellar Object (YSO) jets (Massaglia et 
al. I2005P was carried out with PLUTO, using a simpli- 
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ficd model for the radiative cooling losses, which evolved 
in time only the ionization fraction of hydrogen (c.f., Rossi 
et al. I1997p . This model will be called from now on SNEq 
(Simplified Non-Equilibrium cooling). 

We illustrate a more general treatment of atomic cool- 
ing and evolution of the ionization fraction of the emit- 
ting species, embedded in the PLUTO code as well, for 
use within MHD simulations of astrophysical interest. We 
will call this new cooling function MINEq (Multi-Ion Non- 
Equilibrium cooling). The main advantage of our approach 
is the full ionization state computation during the MHD 
simulation, which allows for better predictions of emission 
line intensities. 

Section [5] contains a general overview of the adopted 
method and implementation of the treatment of radiative 
losses. Then, in Sect. [3J a description of the physics of the 
cooling model can be found, followed in Sect. [J] by the val- 
idations and tests in equilibrium conditions. The numer- 
ical implementation and testing are discussed in Sect. [51 
while in Sect. Owe present some typical astrophysical appli- 
cations. Technical details on ionization-recombination pro- 
cesses and numerical issues are presented in extended form 
in the Appendix. 

2. General overview 

The general characteristics and application ranges of the 
new cooling function added to the PLUTO code are sum- 
marized below. The density limits are those typically en- 
countered in clouds and YSO jets, while the temperature 
range is limited by the highest ionization stage considered 
(at the high end) and the lack of molecular cooling (at the 
low end): 

N e (1(T 2 , 10 5 ) cm" 3 , T G (2 • 10 3 ,2 • 10 5 ) K. (1) 

However, the module is designed to permit later extension 
in terms of applicable parameter range (through adding 
more ion species, or a tabulated cooling function for higher 
temperatures) and physical processes taken into considera- 
tion. 

Flow variables such as density p, velocity v, magnetic 
field B, and total energy E are evolved according to the 
standard MHD equations: 



d(pv) 
dt 



V • (pvv 

OB 
~dt 



T BB T + \ Pt ) 
V x (v x B) 



dE 
~dt 



V-[(E+p t )v-(vB)B] 






(2) 





(3) 





(4) 


Se , 


(5) 



where Se (described later) is a radiative loss term, and pt = 
p+ \B\ 2 /2 denotes the total pressure (thermal + magnetic) 
of the fluid. We assume an ideal equation of state by which 
the total energy density becomes 



P 



\B\ 



r - i 



p- 



(6) 



with r = 5/3 being the specific heats ratio. 

The cooling model accounts for the evolution of 29 ion 
species, namely: Hi, Hn, Hei and Hell, and the first five 



ionization stages of C, N, O, Ne and S. Sulphur, although 
not having an important contribution to cooling, is added 
for diagnostic purposes (line ratios). The ionization network 
employed is larger than in most other MHD codes. 
For each ion, we solve the additional equation 



d{pX^) 
dt 



V • (pX K , 



pS K , 



(7) 



coupled to the original system of conservation laws ([2"|)~ 
(O. In Eq. ([7]) and throughout the following, the first in- 
dex (k) corresponds to the element, while the second in- 
dex (i) corresponds to the ionization stage. Specifically, 
X K ^ = N K i/N K is the ion number fraction, N Ki i is the 
number density of the z-th ion of element k, and N K is the 
element number density. We denote the whole set of ions 
for all possible k and i with X = {X K ,i}. 

The source term S K .i accounts for ionization and recom- 
bination and will be described in the following. The total 
line emission from these species enters in the source term 
Se in Eq. (J5j and should give a good approximation of ra- 
diative cooling for the above conditions (Raga et al. I1997p . 

The system of Eqs. ©-© together with Q is inte- 
grated using the PLUTO code for computational astro- 
physics (Mignone et al. I2007|) . We take advantage of op- 
erator splitting techniques, where the homogeneous part of 
the equations (i.e., with S K) i = Se = 0) is solved sepa- 
rately from the source step. The order of the respective ad- 
vection and source operators (Ti. At and iS At ) is reversed 
every step by keeping the time step At n = At n+1 constant 
for two consecutive integrations to guarantee formal second 
order accuracy. Thus, if U = {p, pv, B, E, pX} is the vec- 
tor of conserved variables, the solution advances from t n to 
t n + At 71 as 



At") = s At "n At " u{t n ) , 



and from V 



At n to t n + 2At n as 



2Af 



n Atn s At " u(t n + At n ) . 



(8) 



(9) 



A new time step, At n+2 , 



is then computed as shown in Sect. 



3. Cooling module description 

We will restrict our attention to the source step only and re- 
mind the interested reader of the original paper by Mignone 
et al. (|2007p for implementation details on the solution of 
the homogeneous MHD equations. 

During the source step, in virtue of operator splitting, 
only internal energy p/(T — 1) and ion fractions X K ^ are 
affected. Density, velocity, and magnetic fields remain con- 
stant with the values provided by the most recent step. 
Thus Eqs. ([5j and (J7J) are treated as a system of ordinary 
differential equations (ODE): 




(r - i)s E 



S K ,i 



(10) 



where k = H, He, C,... labels the element and i = I, II, III,... 
identifies the ionization stage. Equations (jTUJ) must be 
solved for a time increment At n with initial condition pro- 
vided by the output of the previous step (i.e., either an 
advection or source one). 
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Pressure p and temperature T are related by the ideal 
gas equation: 



p = NksT with 



N 



m u fj,(X) ' 



(11) 



where N is the total particle (atoms + electrons) number 
density, ks is the Boltzmann constant, m u is the atomic 
mass unit, and n(X) is the mean molecular weight: 



(12) 



In Eq. (|12j) m K is the atomic mass (in units of m u ) of ele- 
ment k and 7i denotes the number of the ionization stage 
in spectroscopic notation for each ion. B K is the fractional 
abundance of the element. 



3.1. Radiative losses 

Radiative losses are described by the source term Se in the 
energy Eq. (J5J): 



Se = -[ N at N cl A (T, X) + L FF + Li_ 



(13) 



where A(T, X) is the radiative cooling function due to 
collisionally-excited line radiation, Lff denotes the free- free 
(brcmsstrahlung) losses from H + and Hc + , while Zi_r ac- 
counts for the energy lost during ionization/rccombination 
processes. The number densities 7V at and 7V G i are, respec- 
tively, the total atom and electron number densities, readily 
determined from the mass density and the known chemical 
composition of the plasma (by default supposed solar, but 
customizable by the user): 



(14) 



(15) 



Note that iV at does not depend on the ionization state of 
the elements and it should not be confused with N, the 
total particle (atoms + electrons) number density used in 
the equation of state fTTjl. 

Emission lines from Fen, Sin, and Mgil that exist in 
SNEq are added empirically to the energy losses of MINEq 
(without evolving the respective ion species) because of 
their importance at low temperatures. 



3.1.1. Energy loss by collisionally-excited line radiation 

The main contribution to radiative cooling comes from col- 
lisional excitation of low- lying energy levels of common ions, 
such as O and N. In spite of their low abundances, these 
ions make a significant contribution because they have en- 
ergy levels with excitation potentials of the order of kT. 
The total radiative cooling function A(T, X) used in the 
energy source term (Eq. Q3])) is: 



A(T, X) = x ^AN c i, T)B K 



(16) 



where the sums are extended to all ion species and B K 
the fractional abundance of the element n. 



by 



Individual contributions to the different £'s[3 are given 



3 Kj 



(17) 



where Nj is the population of the j-th excitation level; Aji 
are the Einstein A coefficients; and Vji the emission line 
frequency for a transition between levels j and I. We con- 
sider a 5-level atom model to compute the line radiation 
(Ostcrbrock & Fcrland 2005]) by solving for the equilibrium 
populations in each of the excitation levels j = 1 . . . 5: 

N t N cm +^ NiAij = ]T N,N cm +Y^ NjAji , (18) 

W l>3 l& Kj 

which, together with the normalization condition for the 
total number density of the ion, ^2jNj = N K j, can be 

solved for the relative Nj populations in each level. 

The 5-level atom model provides the great majority of 
the emission lines for the considered range of temperatures 
and thus gives a reliable estimation of the total line cooling. 

For most of the ion species, the emission coefficients 
were taken from Pradhan & Zhang (|1999p . The data for 
hydrogen was taken from Giovanardi et al. (|1987p and their 
fit formula. The Cn data comes from Blum & Pradhan 
(|1992p . while for Nil and Nm Chebyshev polynomial fits 
from Stafford et al. (|1994p were used. 

3.1.2. Free-free radiation 

A minor contributor to the cooling rate at moderate tem- 
peratures is the brcmsstrahlung (free-free) radiation, hav- 
ing a continuous spectrum. The rate of cooling in this pro- 
cess by ions of charge Z, integrated over frequency, is ap- 
proximately (Osterbrock & Ferland [2005P 



Lff = 1-42 x 1CT 27 Z 2 T 1/2 N cl N+ 



(19) 



in ergscm _3 s _1 . Because of its abundance, H + dominates 
the free-free cooling, and He + can be included along with 
it since both have the same charge: N + = A^hii + -ZVhcII- 



3.1.3. lonization-recombination losses 

Thermal energy is absorbed by the atom to pass to the 
next ionization stage. During the recombination, a free elec- 
tron is captured and a part of its thermal energy is lost. 
The ionization/rccombination losses are treated similarly 
in MINEq and SNEq, following the method described in 
Rossi et al. PM7| : 



Hie 



-N Q 



(20) 



£l-R = Lj + Lr, 

Lj = 1.27- 10~ 23 y 
L R = 2.39 • 10- 27 VTN HII N cl 

expressed in ergs cm -3 s" 1 . 

The complementary effects of these processes on the 
plasma (the creation/destruction of a free particle) are ac- 
counted through the mean molecular weight, which varies 
together with the total particle number density. 



1 in this section only, k and i will be omitted unless necessary 
to avoid cluttered notation 
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3.2. Ionization network 

Our ionization network can be written in terms of the 
source-term S Kt i mentioned above (Eq. (jTD|) ) : 

S K ,i = N e i X K ^ + ia Ky i+i—X K: i + a K! i)+X Kt i-i£ K! i-i , 

(21) 

where and a K;2 ; are the ionization and recombination 
coefficients of the i-th ion specie of the element k, defined 
as follows: 



N m 



(22) 



^a*?(T), (23) 



where Nun = A^h^hii = -^at-SH^Hii is the number density 
of protons, Nm and iVnei are the number densities of neu- 
tral hydrogen and helium, respectively. The a K ^ and Ck,i 
coefficients arc the transition rates corresponding to the re- 
action mechanisms defined in the Appendix A (note that 
is the total electron-ion recombination coefficient, that 
is dielectronic plus radiative, a cl = a DR + a RR ). 

Since we only consider part of the ions from each ele- 
ment (up to the fourth level of ionization, except for H and 
He), the ionization rate for the highest state will be set to 
zero. This will produce saturation of the ion population in 
this state at very high temperatures, and limit the applica- 
bility of the cooling function. The temperature range can 
however be extended by adding further ionization stages 
for the elements. 

For efficiency purposes, the ionization and recombina- 
tion coefficients on the right-hand side of Eqs. (|2"2")) and 
(|23|) are sampled at discrete values of temperature at the 
beginning of integration and used as lookup tables. 

4. Comparison with equilibrium models 

We perform theoretical line ratios tests to verify the colli- 
sion strengths in the radiative losses. Also, the total cool- 
ing function for an equilibrium ionization balance function 
of temperature (the effective cooling curve) was tested and 
found to be consistent with results obtained with more com- 
plex models. 



4.1. Equilibrium ionization balance 

The equilibrium ionization balance may be used as an ini- 
tial condition for numerical simulations, and also serves for 
testing the ionization/recombination coefficients employed 
in the ionization network. 

The ionization balance for each element at equilibrium is 
computed by setting dX Kyi /dt = (for all ions) in Eq. (fT0)l . 
The equation for the highest ionization stage is replaced by 
the normalization condition, 



K 



(24) 



where K is the highest ionization state taken in considera- 
tion for the element k. Thus, for each element, we solve the 
following system of equations: 

^K,i+l a K,i+l — (Cc,i + a K,i) + -X^.i— lCre,i— 1 = , (25) 



with i — 1, ■ • ■ , K — 1 complemented by Eq. (|24[) . Despite 
its aspect, the previous system of equations is not linear 
since the £ and a coefficients depend on the concentrations 
themselves (see Eqs. (|2"2"j) and (|2"3")l ), so an iterative proce- 
dure must be employed to converge to the correct solution. 

In the particular cases of hydrogen and helium, because 
of the charge-transfer reactions they are involved in, an 
exact treatment would also force £ and a to depend on the 
number densities of all other ions that take part in these 
processes. Considering the very limited influence of these 
reactions on the hydrogen and helium ionization balance, 
we chose to neglect such influences. 

Given an initial guess on the ionization state of the 
plasma, the systems of equations for equilibrium are solved, 
providing new values of iV e i, Nm, and Nu e i- The process 
is repeated until the differences between the old and the 
new solutions are below a certain threshold. The conver- 
gence is rapidly achieved, generally less than five iterations 
are needed for a 10~ 4 — relative threshold. This is ac- 
ceptable, considering that this equilibrium computation is 
typically done on the whole computation grid only once, 
in the beginning of the simulation. In Fig. [TJ we show the 
equilibrium ionization balance as a function of temperature 
for three selected elements. Our results favourably com- 
pare to those obtained by previous investigators-such as 
Sutherland & Dopita (|1993p -with the ionization fractions 
being within 5 — 10% at the same temperature. 

4.2. Line ratios tests 

These tests are useful to verify the emission lines data and 
the level population computation routine (in our 5-level 
atom model). An example is presented here. 

A popular way of estimating the temperatures in 
gaseous nebulae is to use the ratio of spectral line inten- 
sities, such as the lines of Om: 

e(A5007) + e(A4959) e( 1 D 2 -> 3 P 2 ) + e( 1 D 2 -> 3 Pi) 



e(A4363) 



efSo -> 1 D 2 ) 



(26) 

Inserting numerical values of the collision strengths and 
transition probabilities (Osterbrock & Ferland [2005) , the 
ratio becomes: 



R = 



e(A5007) + e(A4959) 
e(A4363) 



8.32 exp 



3.29 x 10 z 
T 



1 + 4.5 x 10" 



r i/2 



, (27) 



for temperatures around 10000K. 

Line ratios computed with the previous formula and the 
results of the 5-level atom model were compared, the dif- 
ferences of less than « 6% being due to the fact that our 
code uses temperature-dependent collision strengths, while 
in the formula above they are assumed constant. 

4.3. Effective cooling 

In Fig. we plot the effective cooling function (in 
erg cm 3 s -1 ) using the ionization fractions computed at 
equilibrium. For the sake of comparison, we also show the 
results obtained with the SNEq model described in Rossi 
ct al. (|1997|) and the Cloudy atomic code, which has a 
large chemical network (see Ferland et al. I1998j) , for similar 
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Fig. 1. Ionization fractions at equilibrium for the five ion- 
ization levels considered for C (top panel), N (middle 
panel), and O (bottom panel). 



plasma conditions. Solar abundances have been assumed 
for all cooling functions, except for the Cloudy Z = 0.3 
case (where the metallicity is only 0.3 times the solar one). 
The SNEq model consists of the emission of 17 most impor- 
tant lines, plus the two-photon continuum and the radia- 
tive losses from ionization/recombination processes. In this 
model, however, only the ionization of H is evolved with the 
integration, the rest of the ions abundances being fixed or 
locked by charge-transfer processes to the ionization state 
of H. 

The results show a good agreement between the newly- 
developed cooling model (MINEq) and the computations 
carried out with the Cloudy code. Chemical composition 
(more extended in Cloudy that is an atomic code) and 
physical processes considered account for the differences. 
The MINEq effective cooling, considering only few met- 
als, generally lies in between the results obtained with the 
Cloudy code for Z = 0.3Zq and Z = Z@. The faster in- 
crease in the peak at 17000K due to hydrogen Lya pre- 
sented by MINEq is due to the different sources of the ion- 
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Fig. 2. Effective cooling curves in the temperature range 
from 10 3 to 2 • 10 5 K, comparison between the results ob- 
tained with MINEq, SNEq, and Cloudy. 



ization/rccombination and emission coefficients (collision 
strengths). 

It can be inferred that while MINEq accounts with good 
accuracy for the cooling losses up to 2 • 10 5 K, SNEq can- 
not follow them above 3 • 10 4 K because it lacks higher 
ionization stages for the atoms. Furthermore, the effective 
cooling obtained with the MINEq model closely reproduces 
the early work of Dalgarno & McCray (|1972p in the tem- 
perature range considered. 



5. Numerical implementation 

In the source step, we advance the system of ordinary dif- 
ferential equations (ODE) given by Eq. ([TO]) in each compu- 
tational zone. For ease of notations, we rewrite the system 
as 

§ -m. m 

where y = {p, X K ^} and f = {(T - 1)Se, S k a} are, respec- 
tively, the vector of unknowns and right-hand sides for all 
possible values of k and i in a given computational cell. 
According to the notations introduced in Eq. ([8]) and ([9]), 
we write the formal solution to (|28|) for a time increment 
At n as y* = S At y°, where the initial condition y° is given 
by the output of the previous step. 



5.1. Integration Strategy 

Accurate numerical integrations of Eq. (f28|) should be car- 
ried out consistently with the different timescales that may 
concurrently co-exist, according to the initial density, tem- 
perature and chemical concentrations. In addition, the sys- 
tem evolution dictated by the local ionization, recombi- 
nation, and cooling rates may proceed considerably faster 
than the typical time scale imposed during the advection 
step. Under some circumstances, this contrast may lead to 
a stiff system of ODE. A common occurrence takes place, 
for instance, when a strong shock propagates in a cold neu- 
tral medium: as the front advances from one computational 
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cell to the next, the ion populations will try to re-adjust to 
the sudden increase in temperature at a rate given by the 
ionization coefficients. At high temperatures (T > 2 x 10 4 ), 
this process may proceed more and more rapidly. 

Nevertheless, these kinds of events are extremely local- 
ized in space since most regions ahead of and far behind 
the shock wave are either close to equilibrium or evolve on 
much slower recombination scales. This suggests some form 
of selective integration by which regions of the flow under- 
going very rapid changes should be promptly detected and 
treated accordingly. We achieve this by 1) detecting poten- 
tial "stiffness" due to large ionization and recombination 
coefficients given by Eqs. (|2"2"]) - (j2"3")) and 2) monitoring, in 
each computational cell, the accuracy through an estimate 
of the local truncation error. We now describe in detail the 
numerical implementation of a dynamically-adaptive inte- 
gration strategy, also shown in Fig. [3] 

At the beginning of integration, we tag a computational 
cell as "non-stiff" if the integration time step satisfies 



At n < 



1 



N e \ max(|C«,i + a Kii \) 



(29) 



where 7V c i is computed in the considered cell. If the pre- 
vious condition holdfQ, we solve Eq. (|28p using an explicit 
method with adaptive stepsize control. Embedded Runge- 
Kutta (RK) pairs simultaneously giving solutions of order 
m and m — 1 arc preferred, since they provide an efficient 
error estimate. The most simple (2,1) pair (m = 2), for ex- 
ample, may be obtained using a simple combination of two 
right-hand side evaluations yielding, respectively, 1 st - and 
2 nd -order accurate solutions y\ and y^: 



y 



y 



° + At n f(y°) 



y 2 = y° + A<'\f (y + ^-/°) . 



At" 



(30) 
(31) 



The difference between the two solutions, y 1 and y 2 , 
estimates the truncation error of the lower order method, 
0{At 2 ) for m = 2. The solution given by y 2 (RK2) is ac- 
cepted only if the error falls below some predefined toler- 
ance e to i (typically 10~ 5 ): 



max 



P --l 

2 



max -X* 



< e to i 



(32) 



where (k, i) extends to all ion species. A more accurate 
Rungc-Kutta (3, 2) pair may be used instead. The condi- 
tion (|52"j) is usually satisfied in regions close to equilibrium 
ionization balance. If Eq. (f3"2"| is not fulfilled we switch to 
an explicit Rungc-Kutta method of order 5 with an embed- 
ded 4 th order solution with coefficients given by Cash-Karp, 
see Press et al. (I1992p . from now on CK45. The adaptive 
strategy provides a 5 th -order accurate solution and allows 
us to split (if required) the full time step At n into a number 
of smaller sub-steps until the condition (f3"2"| is fulfilled in 
each one of them. 

On the other hand, if Eq. (|29|) is not met, explicit time 
marching may potentially become unstable. In such situa- 
tions, integration is carried using a 4 th order semi-implicit 
Rosenbrock method with a 3 rd order embedded error es- 
timation (Ros34 henceforth). Rosenbrock schemes can be 



2 This is, in fact, half the stability limit for the 1 st order 
explicit Euler method. 



considered linearly implicit generalizations of Rungc-Kutta 
methods, the prototype of which is the semi-implicit back- 
ward Euler method. 



(I - At^) ■ (y 1 - y°) = At"/ (y°) 



(33) 



These methods retain stability for large time steps at the 
additional cost of computing the Jacobian matrix J = 
df/dy of the system and performing matrix inversions by 
LU decomposition. Both features are notoriously time con- 
suming for moderately large systems of equations, such as 
the one we deal with here. In Appendix ([Bj we show how 
the Jacobian can be computed using combined analytical 
and numerical differentiation. The full integration strategy 
is schematically illustrated in Fig. [3] 




Fig. 3. Flowchart of the dynamically switching integration 
algorithm for cooling: Runge-Kutta 12, Rosenbrock 34, and 
Cash-Karp 45. Stiffness is detected according to Eq. (f29l) . 



Alternatively, we also found satisfactory results by di- 
viding the whole step At n into smaller ones and by sub- 
cycling with the explicit CK45 scheme. The sub-time step- 
ping strategy has proved to handle the moderate stiffness 
arising at high (T > 10 5 ) temperatures, when the reaction 
rates become large. This makes, in our experience and for 
the tests presented in this work, the explicit scheme com- 
petitive with the semi-implicit method, inasmuch stiffness 
is spatially confined to a small fraction of the computational 
domain. 

Once acceptable solutions y* have been produced in 
every computational zone, we estimate the next time step 
according to the CFL stability restriction and the maximum 
fractional change produced during the radiation step: 



At 



n+2 



min(At adv ,e max At r 



ad J 



(34) 



where, consistently with Eqs. ([8]) and (J9]), the minimum is 
taken over two consecutive time steps and At a( j v is com- 
puted from the CFL condition. The quantity e max spec- 
ifies the maximum fractional change tolerance (typically 
0.01 < e max < 0.1) allowed during the source step. The 
radiative time step At rac i is computed as 



At r 



At n 



max 

xyz 



- 1 



, max 

K,i 



\ V* V® I 

J K,i K,i\ 



(35) 



Note that small values of e max will result in a better cou- 
pling between the advection and source steps, at the cost 
of reducing the overall time step. 

In terms of right-hand side evaluations, the compu- 
tational overheads introduced by the selected algorithms 
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(RK2:CK45:Ros34) arc in the ratio 2 : 6 : 3 in each cell 
for a given time step. The semi-implicit method Ros34 re- 
quires, however, 2 additional right-hand side calls to form 
the Jacobian (see Appcndix[Bj and the inversion of a matrix 
by LU decomposition. This makes Ros34 the most compu- 
tational expensive scheme of integration. 

Nevertheless, extensive testing confirms that only a very 
small fraction of computational zones (usually < 1%) does 
actually require this special, but nonetheless crucial, treat- 
ment. The remaining vast majority of cells can be accu- 
rately evolved using a second order method. On the other 
hand, unconditional use of the CK45 or Ros34 algorithms 
throughout the whole grid leads to a noticeable loss of com- 
putational efficiency with no gain on the overall accuracy. 

5.2. Accuracy comparison 

In order to test the accuracy and efficiency of the selected 
time marching schemes adopted during the source step, we 
consider the evolution of a single parcel of gas departing 
from initial conditions far from equilibrium. This situation 
is typically encountered, for example, when a strong shock 
propagates into a cold medium. Neutral atoms crossing the 
front will suddenly feel the sharp increase in temperature 
and will try to readjust to the new conditions. The ioniza- 
tion timescale will be, most likely, much shorter than the 
typical advection scale. One is interested in performing the 
simulations at the timestep given by Eq. (|34|) . but this can 
violate the condition expressed by Eq. (|29|) . 

We consider two cases in which a single computational 
zone is being evolved in time. Initial parameters have been 
found by running a full shock simulation like the one pre- 
sented in Sect. 16.11 and selecting the computational zone 
showing the most extreme stiff conditions, according to Eq. 
(|29p . We perform a number of time integrations at con- 
stant step size using the Euler, RK2, CK45, and Ros34 
algorithms previously described. Errors are computed with 
respect to a reference solution obtained with the CK45 in- 
tegrator with a very stringent tolerance (10 -8 ) and a small 
timestep (~ 10 -4 of the cooling timescale): 



E»,j \ X K 



X 



rof I 



"V^rci 



(36) 



The errors in temperature are lower in all cases because the 
equations of the chemical network can be, and usually are, 
more prone to very rapid variations (stiffness). 

In the first case (top panel of Fig. QJ, the initial tem- 
perature is set to T = 132 000 K, the initial neutral hydro- 
gen fraction is 22% and the rest of the elements are in the 
highest ionization stage. Under these conditions, the ioniza- 
tion/recombination timescale is r « 2 • 10 3 s, typically much 
smaller than the scale on which hydrodynamical variables 
are transported, At. 

At a fixed timestep At = 5r (sec middle panel in Fig. 
RK2 is less accurate than Euler, this being a typical sign 
of stiffness (Ekcland ct al. 11998(1 . The integrator Ros34 
yields the best accuracy, immediately followed by CK45. 
As the time step is further increased (At = 50t, see bot- 
tom panel in Fig. |4]), CK45 progressively loses accuracy, 
whereas Ros34 keeps the smallest errors. In this case, the 
use of a semi-implicit method clearly reveal its advantages. 

In the second test, we consider a fully-ionized gas (ex- 
cept for hydrogen, AThi « 69%) at low temperature T = 



1.4 
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Fig. 4. Top panel: Temperature and ionization fraction of 
hydrogen evolution; Middle panel: relative errors in ioniza- 
tion fractions, for a fixed timestep At = 5r; Bottom panel: 
same as Middle, for a timestep At = 50r. The error plots 
begin at a time At from the beginning of the simulation. 
Linear time axis up to 10 4 s, logarithmic after. 



10 4 K. For this choice of parameters, the recombination 
timescale is even smaller than before, r sa 10 3 s. Figure [5] 
shows the errors computed with selected integration al- 
gorithms when At = IOOt. The resulting accuracies con- 
firm the trend found for the previous case: low-order, non- 
adaptive time marching schemes are not suitable in con- 
ditions far from equilibrium. On the contrary, Ros34 be- 
ing an adaptive semi-implicit scheme, does not suffer from 
this loss in accuracy and turns out to be the best integra- 
tion method. Explicit adaptive algorithms such as CK45, 
still exhibits somewhat better results than the lower order 
methods. It should also be mentioned that the accuracy 
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of explicit schemes may be further improved if step sub- 
cycling is used. 




0.1 



Fig. 5. Top panel: Temperature and ionization fraction of 
hydrogen evolution; Bottom panel: relative errors in ioniza- 
tion fractions, for a fixed timestep At = IOOt. Linear time 



axis up to 10 s, logarithmic after. 



We conclude that, for the slow varying regions of the 
MHD simulation, RK2 (or RK3) can be a good choice, 
while for very fast varying regions a higher-order integrator 
with time step adaptivity (as CK45) or even an implicit one 
(Ros34) become necessary. Also, large temperatures are not 
a necessary condition for the system of equations to become 
stiff, since this can also happen at relatively low tempera- 
tures when the ionization/recombination rates are high. 



6.1. Propagating shocks 

It is interesting to see the difference radiative losses make 
in the dynamical evolution of a propagating shock. A first 
scries of tests were made in ID, with an initial perturba- 
tion in pressure, density, and velocity that propagates in a 



stratified medium of T, 



pre 



1 00QK and becomes a shock. 



The pre-shock density in the external medium is 

/ \ x o 

P0[X) = PO-T- 2 ! 



(37) 



where x is the spatial coordinate and the departure density 
Po corresponds to a particle number density A*o = 10 5 cm~ 3 . 
This density distribution should approximate well the den- 
sity in an expanding jet. The initial perturbation is set in 
such a way that only one shock forms instead of the usual 
pair of forward/reverse shocks. The setup is described in 
detail in Massaglia ct al. (2005). The ID simulation was 
run on a domain of length L = 4 x 10 15 cm, with a resolu- 
tion of 1.4 x 10 11 cm, and the initial velocity perturbation 
had an amplitude Av = 30kms _1 . 

In Fig. [6l a comparison is made between the evolution 
of the formed shock in the absence of cooling, with SNEq 
and the evolution with MINEq. Plots of density and tem- 
perature are presented at three evolutionary instants in the 
propagation. As it results from the plots in Fig. [51 the shock 
dynamics arc heavily influenced by radiative cooling. The 
shock propagation velocity decreases almost twice in the 
simulations with cooling with respect to the adiabatic ones. 
The smooth decrease in temperature after the shock front 
in the adiabatic simulation is replaced by a much sharper 
one in the simulations with cooling. 




Fig. 6. Density logarithmic profiles (top row) and temper- 
ature profiles (bottom row), at three evolutionary stages: 
adiabatic (dotted line), with SNEq (dashed line), and 
MINEq (solid line) cooling. 



6. Astrophysical applications 

We now apply the newly-developed cooling function to 
problems of astrophysical interest. First, we consider a sin- 
gle, one-dimensional radiative shock propagating in a strat- 
ified medium with decreasing density. Then, an example of 
application of the first setup for the computation of emis- 
sion line ratios is shown. Finally, a study of the dynami- 
cal evolution of a jet with varying ejection velocity in two- 
dimensional axial symmetry is presented. 



The differences between the density plots obtained with 
the two cooling models arc quantitatively moderate, while 
the shapes are similar. The maximum temperatures at- 
tained are very close in the two cooling models. Overall, 
the dynamical differences that appear between the use of 
the two cooling models are small, but the line intensity ra- 
tios are very sensitive to density/temperature conditions so 
the differences may result in moderate amplitudes. 

The test in equivalent configuration was also performed 
in a 2D slab. The results were, as expected, very similar to 
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the ones from the ID simulations, with somewhat smoother 
curves due to the lower resolution employed. It results that 
for simulations of propagating shocks like the one described 
it is very important to include the radiative cooling losses, 
which heavily influence the dynamics. A simplified treat- 
ment of these losses can be sufficient for studies on the dy- 
namics, while for the computation of emission line maps the 
detailed (MINEq) approach is more suitable. An example 
is presented in the following section. 



6.2. Emission lines 

The computation of emission line ratios from numerical 
simulations is of great importance for the field to compare 
to observations and to discriminate between theoretical as- 
trophysical models. 

In a simple ID setup, one of the ways to model a YSO 
jet and to estimate the emission is the following. Supposing 
that the emission comes from shocks inside the jet, the 
propagation of a shock resulting from a velocity fluctuation 
is simulated in the frame of reference of the jet material. 
The emission in the chosen lines is computed and aver- 
aged over a space region corresponding to the resolution of 
the observational data (in our example 10 15 cm) for each 
evolutionary step. Then, the resulting averaged line ratios 
are plotted at space points corresponding to the transport 
speed of the jet material, set to 150fcm • s _1 . 

The computation was done for the setup presented in 
Sect. 16.11 The resulting plot, presented in Fig. [7j has also 
the x axis converted in arcseconds (for a distance D = 
140pc) and represents the emission of a jet in the assump- 
tion that the emission comes from internal shocks formed 
due to initial jet velocity variability. We present ratios be- 
tween the forbidden emission lines of Oi AA6300A + 6364A, 
Nil AA6548A + 6583A and Sn AA6716A + 673lA. Such syn- 
thetic emission line ratios can be directly compared with 
observations of YSO jets. 



l.o : 
0.5 : 




13 3 4 

Distance (arcsec) 

Fig. 7. Line ratios with MINEq cooling, for the propagating 
shock described in the previous section. 



The main advantage in using MINEq for creating syn- 
thetic emission maps is that non-equilibrium ionization bal- 
ance for the elements are provided. The computation of the 
emission in selected lines is done in post-processing, with 
routines distributed together with the code. 



6.3. Jet propagation 

Observations of YSO jets that show series of emission knots 
along their length pointed out that simple steady-state 
models cannot explain their morphology. These knots have 
been interpreted in the literature as due either to the non- 
linear evolution of Kelvin- Hclmholtz instabilities set at the 
jet-ambient interface (Micono et al. I2000|) or to velocity 
variabilities of the beam (Internal Working Surfaces, see for 
example Raga et al. I2002|) that, during their propagation, 
steepen into shocks. The latter scenario has been chosen as 
a possible astrophysical application of the cooling module. 

In the present case, we consider a variable jet in 2D 
cylindrical geometry propagating into a uniform ambient 
medium with particle number density n a = 200 cm~ 3 and 
temperature T a = 2 500 K. The beam is injected at the 
z = 0, r < Rj (i?j = 2.5 x 10 15 cm) boundary with higher 
density (nj = 5 n a ) than the background. The mean jet 
velocity is Vj = 110 km s^ 1 with sinusoidal oscillations of 
amplitude Av = 25 km s™ 1 and a period of r = 50 yrs 

A purely toroidal magnetic field is injected at z = 
along with the beam, following the simple configuration 
described in Lind & al. (|1989|) : 



B m - for 0<r<a, 
a 



B <t>( r ) = { B m - for a<r<R, 



(38) 







otherwise , 



where B m and a are the magnetization strength and ra- 
dius. Demanding pressure equilibrium at the jet inlet, 
d(p + B^)/dr = —B 2 /r, one recovers the pressure profile 
inside the beam (r < i?j), 



P(r) =Po~ B 2 n min ( 1, 



(39) 



where po is corresponds to a central temperature Tq = 
10 000 K. Finally, the magnetization strength, B m is pre- 
scribed from the plasma ft parameter, defined in terms of 
the averaged beam pressure: 



ft 



2 



/(f 3 p{ r ) r dr a 2 
f* s rdr 



4 



,_Pp_ 

' B ' 2 



(40) 



from which one can easily recover B m . For the present com- 
putation, we set a = 0.6i?j and [3 = 1. This choice of pa- 
rameters is similar to the ones found by Masciardi & Raga 
(|200ip in their attempt to model the curved HH 505 jet. 

Numerical integration is carried out with the PPM 
method and the HLLC Riemann solver of Li (|2005p . We 
use 30 zones on the jet radius and the domain extends, in 
beam radii, from to 10 in the radial direction and from 
to 60 in the longitudinal direction. Free outflow is assumed 
across the outer boundaries, whereas reflecting boundary 
conditions hold at the axis (R = 0) and outside the jet in- 
let (z = 0, r > Rj). A smoothing function is introduced for 
all variables at the transition between the jet material and 
the external medium to avoid the formation of an unphys- 
ical high temperature low-density layer around the jet. 

We perform a set of three simulations, by adopting i) 
a tabulated cooling function, ii) simplified treatment of ra- 
diative cooling losses (SNEq), and iii) the detailed cooling 
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Fig. 8. Tabulated Cloudy (top panel), radiative with SNEq cooling (middle panel) and MINEq cooling (bottom panel) 
for the jet simulation. In each panel, the upper and lower halves show density and temperature in log scale, respectively. 
The different shades are normalized between p m i n and /9 max (for density), T ro j n and T max (for temperature). 



treatment. The tabulated cooling function simply consists 
in adding a source term to the energy equation, given by 
the Cloudy Z=l cooling curve (presented in Sect. I4.3j) as 
a tabulated function of temperature. This cooling imple- 
mentation does not follow the ionization balance of any 
element, but as a standard procedure for this kind of ap- 
proach, the effective cooling function is multiplied by the 
particle density squared to obtain energy losses per unit 
volume. Results at t w 500yrs are shown in Fig. JHJ). 

The pulsed initial jet velocity produces, as expected, a 
number of intermediate shocks propagating along the jet 
with typical post-shock temperatures in the range 15 000 
25 000 K. The morphology is similar between the SNEq and 
MINEq runs, since at these temperatures the two cooling 
losses are comparable. Larger deviations are observed close 
to the head of the jet, where temperatures are higher and 
the two cooling functions exhibit larger differences. Overall, 
while SNEq and MINEq give similar results, radiated losses 
are higher for the tabulated Cloudy cooling, as can be in- 
ferred by the reduced lateral expansion of the cocoon (this 
can be expected considering the effective cooling curves in 
Fig.©. 

Nevertheless, temperatures at the jet head are highest 
for the tabulated cooling case. In order to understand this 



apparently unexpected result, we have performed a system- 
atic comparison between the Tabulated and SNEq cool- 
ing functions by means of supplementary simulations (not 
shown here). Our results demonstrate that at high tem- 
peratures (> 4 • 10 4 K) and low ionization, the SNEq line 
emission becomes larger than the equilibrium Cloudy one. 
In this case, it is crucial that the SNEq line emissions de- 
pend on the electron number density and that this density 
is dynamically computed from the non-equilibrium ioniza- 
tion of H. In the tabulated case, an equilibrium ionization 
balance is implicitly assumed. This confirms that the max- 
imum temperatures can be an indication of the maximum 
cooling losses attained locally, but not on the overall cooling 
losses. It must also be noticed that the resolution of these 
simulations is still low to resolve the post-shock zone at the 
jet head, so the maximum temperatures observed may be 
subject to large uncertainties. 

We conclude that even a simple cooling like SNEq, 
evolving only the hydrogen fraction, is a much better ap- 
proximation than using tabulated cooling losses. 

The ionization fractions computation is very important 
when it comes to producing synthetic emission maps in 
various emission lines to be compared with observations. 
In Fig. ©, the fractions of Nil and On are presented, dy- 
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namically computed by MINEq and alternatively computed 
from SNEq considering them fixed by the hydrogen ion- 
ization through charge-transfer (see Osterbrock & Ferland 
2005). The differences are moderate and can result in vari- 
ations of 20 — 30% in the emission lines computation. 

The steep gradients and transition regions forming im- 
mediately behind the shock front are crucial in determin- 
ing the emission properties, e.g. line intensity ratios. For 
this reason they need to be accurately resolved (Massaglia 
et al. I2005|) . However, at the resolution employed here 
(300 x 1800), only the general physical evolution can be 
captured. Considerably higher resolution is required and 
this can be efficiently achieved only through adaptive mesh 
refinement simulations. This issue will be the subject of 
forthcoming works. 

7. Conclusions 

After a series of tests and validation, the detailed treatment 
of radiative losses MINEq is now implemented in the MHD 
code PLUTO. The choice of the integration technique and 
the optimizations are, as far as we know, unique at the 
time of this writing and provide a high degree of accuracy 
in ion species abundances and radiative losses computation. 
Both theoretical and technical aspects of the current imple- 
mentation, as well as testing process and applications were 
discussed in the previous sections. 

A major advantage obtained by using MINEq, com- 
pared to previously employed SNEq, is that the line emis- 
sion can be computed in conditions of non-equilibrium ion- 
ization for all species, more likely to be encountered in sit- 
uations of rapid changes, as it is the case of shock waves. 

As shown by the tests presented, the choice of the cool- 
ing model between MINEq and SNEq, has an increasing 
effect on the structure formation whenever temperatures 
exceed 25 000 K, and has an important influence for the 
ionization fractions that reflects in the emission line com- 
putations. For a preliminary dynamical study, a tabulated 
cooling function that does not integrate any ion specie can 
be employed, but the relatively low computational cost of 
a cooling that evolves the hydrogen ionization in a time- 
dependent fashion makes the latter advisable in most cases. 
Whenever the purpose is the computation of emission line 
ratios, employing a detailed cooling like MINEq produces 
more reliable results as the ionization fractions are followed 
in non-equilibrium conditions, twhich arc likely to be en- 
countered in the real astrophysical objects. 

An important feature of the cooling model implemen- 
tation is that it is upgradeable, more ion species and other 
processes can be added (increasing also the temperature 
range of application). Also, it can be used as starting point 
for the integration of atomic chemistry and cooling pro- 
cesses to other MHD codes. 

The newly-developed cooling function provides a pow- 
erful tool for investigating the stellar jets and gaseous neb- 
ulae. It is, in the current configuration, suited for the study 
of radiative shocks in stellar jets. High resolution, adaptive 
numerical simulations to predict line emission in YSO jets 
will be subject of forthcoming works. 

Appendix A: lonization/recombination processes 

The following processes are taken into consideration: colli- 
sional ionization, radiative, and diclcctronic recombination, 



charge-transfer with H and He. These processes enter the 
ionization/recombination coefficients defined and used in 
Sect. H 

A.l. Collisional ionization 



We use the Voronov (|1997|) data to estimate the collisional 
ionization rates with the analytical formula: 



Tip. Tjl/2 

■■ x + u ■ £/K - e " U ' 



(A.l) 



where U = AE/T and A, P, AE, X, and K parameters are 
listed in Table 1 of the cited paper. T and AE are expressed 
in eV, and £ co11 in cmV 1 . 

The actual number of ionizations in unit time and unit 
volume will be: 



dt 



(A.2) 



where Ni is the total number density of atoms in the lower- 
ionization state, and N e \ the electron number density. 

A.2. Radiative recombination 

The total radiative recombination rates are taken from 
Pcquignot et al. (199T|). The total recombination rate is 
fitted with the analytical formula: 



a RR = 10 -13 z 



at" 



1 + ct° 



(A.3) 



where z is the ionic charge and t = l0~ 4 T(K)/z 2 . The four 
parameters a, b, c, and d are given in Table 1 in the cited 
paper. 

The resulting cv RR is expressed in units of cm 3 s _1 . 

A.3. Dielectronic recombination 

The diclcctronic recombination process proceeds as 



.4 



+111+1 



Al 



Al 



hv. 



(A.4) 



where p stands for a state of the m+1 times ionized element 
A, and a and b represent an auto-ionizing and a true bound 
state of the next ionization stage. 

From Nussbaumer & Storey Q1983p . the dielectronic re- 
combination rates are fitted by the analytical formula: 



„DR 



10 



-12 



[j+b + ct + ct 2 )t- 3 / 2 exp(^j, (A.5) 



where t — T(K)/1Q A K and cv DR is expressed in cm 3 s _1 . 
The coefficients are given in a table from the cited paper. 
We used the data from Nussbaumer & Storey (|1983p for 
the C, N, and O ions. 

A.4. Total electron - ion recombination 



For the He, Nc, and S ions, data from Kato & Asano (|1999[) 
was used for the total recombination coefficient (radiative 
+ dielectronic). These arc tabulated values that we inter- 
polate in our temperature range. 
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Fig. 9. Fractions of N II (top panel) and O II (bottom panel) for the jet simulation. In each panel, the upper and lower 
halves show the results with SNEq and MINEq, respectively. The different shades are normalized between the minimum 
and maximum values of the fractions. 



A. 5. Charge transfer with H 



Appendix B: Jacobian Matrix 



The charge transfer (exchange) reactions with H are reac- 
tions of the form: 



A+ n + A +( - n - 



SE 



(A.6) 



The direct reaction is called charge-transfer recombination 
and the inverse charge-transfer ionization (C HI1 )- 

We took the data for charge transfer with H from 
Kingdon & Ferland (|T996). The recombination/ionization 
rate is fitted by the analytical formula 



HI aHII 



C mi = <[i 



(A.7) 



where t± = T(K)/10 4 K and the parameters a, b, c, and d 
are listed in Tables 1 and 3 from the cited paper. 

A.6. Charge transfer with He 

The charge transfer (exchange) reactions with He are reac- 
tions of the form: 



A +n + He -> A+( n " 



He 



(A., 



The solution of implicitly linearized equations such as Eq. 
(|33[) or the Rosenbrock scheme requires the expression of 
the Jacobian matrix of the system of equations given by 
(fTU|) . Using the definition, 



J 



= df(y) 
dy 



( dX dx\ 

dX ~dp~ 

dp dp 

\dX dp~ ) 



(B.l) 



where y = {p, X} and f(y) = y. Partial derivatives are 
computed using combined analytical and numerical differ- 
entiation. We note in the first place that, for practical rea- 
sons, the right-hand side f(p,X) is better expressed in 
terms of temperature and ionization fractions, that is 



(B.2) 



f[p,Xj=g[T,X 

where T, p, and X are related through 

pm u n(X) 



T 



(B.3) 



We took the data for charge transfer with He from Wang This allows us to compute partial derivatives with respect 
et al. (i200ip and references herein. The recombination rate to the ion fractions using the chain rule, 



is fitted by the analytical formula 

a HcI = at\[\ + cexp {dt 4 )} : 





df 




dg 


H 


dg 


8T 


(A.9) 


dX^^ m 






T 


' 0T 


x dX^ m 



(B.4) 



where £4 = T/10 4 K and the parameters a, b, c, and d are 
listed in tables available on-line. 



where, using Eq. (|B.3[) and the definitions of the mean 
molecular weight, Eq. (|12|). we can express the second term 
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on the right as 
dT 



dg 

dT 



dj_ 

dp 



P 

x 



1 dfi 



N 



(i d(i L 



\x D dX^ 

(id 9Xc, 



(B.5) 

where the term is square brackets is simply d(x(X)/dX^ m 
whereas (in and (Id are, respectively, the numerator and 
the denominator of the mean molecular weight. 

The explicit dependence on T and X in our ionization 
network, Eq. (|2ip . is made clear by rearranging terms as 



X K .i — L K ,i (T, X) X Kt i—i — C K .i (T, X) X Kt i+ 

+ Rn,i (T, X) X Kt i+l , 



(B. 



for each element's ions. The coefficients L Kj i, C Kt i, and i? Kj i 
are expressed by sums of functions depending on either T 
or X: 



L Kl i — L a Kl N c \ + L Ki X m + L c Ki Xn c i + L d Ki 
C K ,i = C%iN el + C b t X m + C^ iX He i + 0% , 



R 



d 

K*i 



(B.7) 
(B.8) 
(B.9) 



where the L' Ri 's, C'^'s, and R'^'i's depend on T only 
whereas JV e i, given by Eq. (JT5J) , depends on X only. 
Focusing on the Jacobian sub-matrix dX/dX, we evalu- 
ate the first term in Eq. (|B.4|) as 



dX R>i 



dX, 



dR K .j 

dX^ m X K ^i 



9X^^ m 

dC Kti 
dX, 



X K ,i~i- 

x K ,,+ 



(B.10) 



where 8i^ m is the Kroncckcr delta symbol and 
dL K i 



OX, 



dC K , 



dX, 



dR Kti 
dX, 



K, l {T)N lm A i +L b ^{T)5 H1 , im + 
+L c ^ i {T)8 H ei^ m 

R a K ,i(T)N lm A 6 + & Kti (T)5 HI ,t m + 
+R c K ,i( T ^ HeI ^ m 



(B.ll) 



(B.12) 



(B.13) 



In the previous equations we made use of the fact that 
dN el /dX K)i = N( li -l)B K . 

The last row of the Jacobian involves derivatives of the 
cooling function with respect to X: 



dp 



dX^^ m 

where 
<9A 



dh 



FF 



dX, 



(B.14) 



dC 

C^ m B s + Nj m B s ]T .V, .,— ^/»", , (B.15) 



dX, 



dX, 



dC K .i/ dX^^m is found numerically and the remaining terms 
are found by straightforward differentiation of the energy 
losses due to ionization-rccombination and brcmsstrahhmg. 

Finally, partial derivatives with respect to pressure 
needed in Eq. (|B.5[) and in the last column of J are com- 
puted numerically using a centered approximation: 



(B.16) 



df a f(p{l + e),X)-f(p(l-e),X) 
dp ep 

where e is a small parameter (typically e = 10 ). 
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